library(datasets)
data(iris)
Avant de débuter des analyses, il est important de se familisariser avec son jeu de données afin d’avoir une idée de sa structure. Cette étape permet d’identifier des motifs, des tendances ou des relations préalablement aux tests statistiques, puis de vérifier la qualité des données. Elle permet aussi d’évaluer si les données ont besoin d’être transformées avant de procéder aux analyses.
Pour vérifier le nombre de lignes et de colonnes de votre dataset:
dim(iris)
## [1] 150 5
Pour voir les 6 premières lignes du jeu de données:
head(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
Pour connaître le nom des colonnes:
names(iris)
## [1] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width" "Species"
Pour vérifier la classe de chaque variable dans le dataset:
str(iris)
## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
Pour obtenir certaines statistiques descriptives de base, telles que le minimum, le maximum puis la moyenne (variables continues), ainsi que le nombre d’observations (variables catégoriques):
summary(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
## 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
## Median :5.800 Median :3.000 Median :4.350 Median :1.300
## Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
## 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
## Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
##
Pour oconnaître les niveaux d’une variable catégorique:
levels(iris$Species)
## [1] "setosa" "versicolor" "virginica"
La fonction plot() est une fonction versatile qui permet de créer une grande variété de figures avec les données brutes. Dépendemment du type de donnée et des arguments fournis à la fonction, plot() peut générer plusieurs types de figures comme:
Une seule matrice de scatterplots figure peut être créée avec toutes les variables continues.
plot(iris)
Pour créer un scatterplot de la relation entre une variable spécifique et une autre, il suffit d’entrer celles-ci dans la fonction, séparées par le signe ~. Par exemple, si on veut visualiser la relation entre la longueur des sépales et celle des pétales dans le jeu de données iris, on peut l’écrire ainsi:
plot(iris$Sepal.Length~iris$Petal.Length)
À l’aide de la fonction boxplot(), on peut visualiser la dispersion d’une variable selon dles groupes d’une variable catégorique. Ce type de figure permet aussi d’identifier rapidement des valeurs aberrantes ou des anomalies.
boxplot(iris$Petal.Length~iris$Species)
La fonction hist() génère un histogramme de la variable qui y est précisée.
hist(iris$Sepal.Width)
Les figures sont des représentations visuelles des résultats. Elles rendent la lectures des résultats principaux plus facile et permettent de mettre en évidence des tendances ou motifs intéressants. On veut pouvoir la comprendre sans avoir à osciller entre la figure et le texte. Elles doivent contenir:
Un titre descriptif: les variables mesurées, les unités de mesure, le nom commun et en latin du taxon (si nécessaire). Le titre doit fournir assez d’information pour qu’on comprenne le contexte de la figure sans devoir se référer au texte.
Titre des axes: doivent comprendre la variable et son unité de mesure
Barres d’erreurs: les inclure pour indiquer l’écart-type autour de la moyenne
Légende: une petite légende peut être nécessaire pour distinguer les traitements (couleur, type de ligne, etc.)
Indice de significativité: la valeur p d’une relation, des astérisques au-dessus de graphes à barres pour indiquer si la relation est significative, etc.
Les tableaux sont un bon choix pour présenter de l’information numérique détaillée. Ils présentent habituellement des résultats plus complexes qui seraient trop encombrants à inclure dans une figure ou dans le texte. Généralement, si les données ne peuvent être présentées en une ou deux phrases, un tableau est nécessaire. Les lignes et les colonnes doivent contenir le nom de la variable ainsi que l’unité de mesure. Un tableau résumé des statistiques peut inclure, par exemple, la moyenne, l’écart-type, les intervalles de confiances, les degrés de liberté, la valeur p et autres statistiques (comme la F value).
Les légendes servent à compléter l’information qui est présentée dans le tableau ou la figure. On y retrouve entre autre les tailles d’échantillons (n), valeur p, des descriptions des abréviations, la méthode de collection, le nombre de réplicats, etc.
Les données ou les figures qui ne contribuent pas directement à l’histoire principale de votre rapport peuvent être rassemblées en annexes. Dans cette section, on peut retrouver, par exemple, des cartes du site d’échantillonnage, du matériel utilisé pour récolter les données, des tableaux avec davantage d’informations sur les modèles (pas de capture d’écran du summary!), etc.
summary(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
## 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
## Median :5.800 Median :3.000 Median :4.350 Median :1.300
## Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
## 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
## Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
##
La fonction de base plot() peut être utile pour visualiser rapidement un jeu de données, par exemple avec un histogramme. Par contre, vous verrez que son utilisation peut devenir limitée lorsqu’il s’agit de réaliser des figures plus complexes ou simplement modifier certains paramètres graphiques.
C’est pour cela que nous suggérons d’utiliser la fonction ggplot(), du package “ggplot2”. Cette fonction permet de réaliser des graphiques de façon plus intuitive et permet de les mettre en page plus facilement qu’avec la fonction plot(). Elle est aussi mieux documentée, il est donc plus facile de comprendre et utiliser les différents aspects de la fonction. Google est d’ailleurs votre meilleur allié pour la réalisation de vos figures!
Contrairement à la fonction plot(), ggplot() fonctionne par couches. Une figure ggplot() commence avec la fonction ggplot(). Elle sert à “préparer” la figure: on spécifie le jeu de données à utiliser, puis on choisit les variables qui formeront nos axes.
La fonction ggplot() nécessite deux arguments: le dataset (jeu de données), puis l’argument aes(). Ce dernier nous permet d’assigner des variables du dataset aux composantes du graphique (par exemple, les axes x et y). Voici un exemple, toujours avec le jeu de données Iris. Nous allons tester si la largeur des sépales diffère entre les espèces.
(Revoir l’atelier 3 pour l’explication de l’ANOVA, incluant les postulats et critères d’utilisation, que je passe ici)
library("ggplot2")
head(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
boxplot(Sepal.Width ~ Species,
data = iris)
On produit d’abord le modèle linéaire, et le test d’ANOVA, avec Sepal.Width comme variable dépendante et Species comme variable indépendante.
modele.SW<-lm(Sepal.Width~Species,data=iris)
#Test anova:
anova(modele.SW)
## Analysis of Variance Table
##
## Response: Sepal.Width
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 11.345 5.6725 49.16 < 2.2e-16 ***
## Residuals 147 16.962 0.1154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(modele.SW)
##
## Call:
## lm(formula = Sepal.Width ~ Species, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.128 -0.228 0.026 0.226 0.972
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.42800 0.04804 71.359 < 2e-16 ***
## Speciesversicolor -0.65800 0.06794 -9.685 < 2e-16 ***
## Speciesvirginica -0.45400 0.06794 -6.683 4.54e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3397 on 147 degrees of freedom
## Multiple R-squared: 0.4008, Adjusted R-squared: 0.3926
## F-statistic: 49.16 on 2 and 147 DF, p-value: < 2.2e-16
# [[[ Commentaire anova ]]]
Ensuite, on effectue le test post-hoc Tukey:
compSW<-aov(Sepal.Width~Species,data=iris)
TukeyHSD(compSW)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Sepal.Width ~ Species, data = iris)
##
## $Species
## diff lwr upr p adj
## versicolor-setosa -0.658 -0.81885528 -0.4971447 0.0000000
## virginica-setosa -0.454 -0.61485528 -0.2931447 0.0000000
## virginica-versicolor 0.204 0.04314472 0.3648553 0.0087802
summary(compSW)
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 11.35 5.672 49.16 <2e-16 ***
## Residuals 147 16.96 0.115
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# [[[ Commentaire Tukey ]]]
On peut ensuite construire notre figure. D’abord, la fonction ggplot(). On spécifie le jeu de données, puis l’argument aes(), ici les axes x et y:
ggplot(iris, aes(x=Species, y=Sepal.Width))
Où sont les données? Il faut les ajouter!
C’est comme ça que ggplot() fonctionne. On crée la base de notre figure, puis on y ajoute les données à l’aide des fonctions geom_*. Par exemple, geom_points nous permet d’ajouter des données sous forme de scatterplots (donc des points), tandis que geom_line nous permet d’ajouter une ligne. Dans notre cas, nous voulons montrer nos données sous forme de boxplot, donc nous utilisons geom_boxplot. Il ne faut pas oublier d’ajouter un + après chaque fonction:
ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_boxplot()
Pour ajouter/modifier d’autres éléments, il faut ajouter des couches. Ainsi, on peut changer les titres d’axes (labs pour labels), changer le thème de notre graphique (theme_bw est plus minimaliste), retirer la grille et centrer le titre:
ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_boxplot()+
labs( x="Espèce",
y="Largeur des sépales (cm)",
title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
On peut aussi décider de changer les couleurs. Ce n’est pas toujours nécessaire ou pertinent, si ça n’ajoute pas d’information à la figure. Il faut faire attention de ne pas surcharger la figure avec le design esthétique, l’important c’est de focuser sur le message que l’on veut transmettre avec notre figure, et ne pas ajouter d’éléments qui sont distrayants. Mais pour l’exercice et se familiariser avec les paramètres graphiques, essayons:
ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_boxplot(fill = "salmon", color = "red")+
labs( x="Espèce",
y="Largeur des sépales (mm)",
title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
On a donc changé le fill (intérieur) et color (trait) de nos boxplot, à même la fonction boxplot(). On pourrait aussi changer la couleur des boîtes en fonction de l’espèce d’iris (qui est aussi notre variable x). Dans ce cas, il faut ajouter un argument aes() à la fonction boxplot(), pour assigner notre variable à un aesthetic. Ce sera notre variable Species:
#Fill:
ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_boxplot(aes(fill=Species), show.legend = FALSE)+
labs( x="Espèce",
y="Largeur des sépales (mm)",
title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
#Color:
ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_boxplot(aes(color=Species), show.legend = FALSE)+
labs( x="Espèce",
y="Largeur des sépales (mm)",
title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
Il y a aussi des palettes de couleurs prédéfinies dans R, qu’on peut utiliser avec la fonction scale_fill_brewer (j’ai choisi la palette “Accent”):
ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_boxplot(aes(fill=Species), show.legend=FALSE)+
scale_fill_brewer(palette="Accent")+
labs( x="Espèce",
y="Largeur des sépales (mm)",
title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
#Qu'est-ce qui arrive si on change Fill pour Color, dans boxplot()? Il faut aussi changer scale_fill_brewer() pour scale_color_brewer():
ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_boxplot(aes(color=Species), show.legend=FALSE)+
scale_color_brewer(palette="Accent")+
labs( x="Espèce",
y="Largeur des sépales (mm)",
title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
Puisque ggplot() fonctionne par couche, on peut aussi choisir d’ajouter les données brutes à notre figure. Par exemple, on ajoute avant geom_boxplot() une autre couche de données avec geom_point(), dans lequel on peut spécifier, encore une fois, l’argument aes() pour assigner la variable qui sera colorée. L’argument position définit le niveau de “jitter”, c’est-à-dire des points un peu “éparpillés” sur l’axe x. Ensuite, on pourrait ajouter l’argument alpha=0.5 dans la fonction boxplot(), pour les rendre semi-transparents:
ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_point(aes(color=Species), position = position_jitter(0.2), shape = 16, size = 2.5, show.legend=FALSE)+
scale_color_brewer(palette="Accent")+
geom_boxplot(alpha=0.5)+
labs( x="Espèce",
y="Largeur des sépales (mm)",
title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(hjust = 0.5))
On peut aussi ajouter des lettres pour montrer nos résultats de comparaison (voir https://statdoe.com/one-way-anova-and-box-plot-in-r/ pour la méthode). En gros, la fonction multcompLetters4 assigne les lettres en fonction des groupes du test de Tukey, puis on constuit une table qui contient ces lettres (LettresSW dans l’exemple), la variable Species, ainsi que la position de la lettre (en haut du quartile pour chaque boîte/valeur de Species).
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(multcompView)
compSW<-aov(Sepal.Width~Species,data=iris)
tukeySW <- TukeyHSD(compSW)
print(tukeySW)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Sepal.Width ~ Species, data = iris)
##
## $Species
## diff lwr upr p adj
## versicolor-setosa -0.658 -0.81885528 -0.4971447 0.0000000
## virginica-setosa -0.454 -0.61485528 -0.2931447 0.0000000
## virginica-versicolor 0.204 0.04314472 0.3648553 0.0087802
cld <- multcompLetters4(compSW, tukeySW)
LettresSW <- group_by(iris, Species) %>%
summarise(mean=mean(Sepal.Width), quant = quantile(Sepal.Width, probs = 0.75)) %>%
arrange(desc(mean))
cld <- as.data.frame.list(cld$Species)
LettresSW$cld <- cld$Letters
print(LettresSW)
## # A tibble: 3 × 4
## Species mean quant cld
## <fct> <dbl> <dbl> <chr>
## 1 setosa 3.43 3.68 a
## 2 virginica 2.97 3.18 b
## 3 versicolor 2.77 3 c
On ajoute nos lettres au boxplot avec la fonction geom_text(), avec comme argument aes() la table que l’on vient de créer. On peut aussi ajouter la statistique de test avec la fonction annotate(), en utilisant les coordonnées x et y dans notre figure.
ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_point(aes(color=Species), position = position_jitter(0.1), shape = 16, size = 2.5, show.legend=FALSE)+
scale_color_brewer(palette="Accent")+
geom_boxplot(alpha=0.5)+
labs( x="Espèce",
y="Largeur des sépales (mm)",
title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(hjust = 0.5))+
geom_text(data = LettresSW, aes(x = Species, y = quant, label = cld), size = 5, vjust=-1, hjust =-5)+
annotate("text", x = 3, y = 4.5, label = "ANOVA, F(2)=49.16, p<0.001")
Pour bien faire, on pourrait même changer les étiquettes en x. Avec la fonction scale_x_discrete, j’ajoute “I.” pour iris à chaque étiquette et je les mets en italique, puisque ce sont des noms latins, avec l’argument axis.text.x dans la fonction theme().
ggplot(iris, aes(x=Species, y=Sepal.Width)) +
geom_point(aes(color=Species), position = position_jitter(0.1), shape = 16, size = 2.5, show.legend=FALSE)+
scale_color_brewer(palette="Accent")+
geom_boxplot(alpha=0.5)+
labs( x="Espèce",
y="Largeur des sépales (mm)",
title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(hjust = 0.5), axis.text.x = element_text(face="italic"))+
geom_text(data = LettresSW, aes(x = Species, y = quant, label = cld), size = 5, vjust=-1, hjust =-5)+
annotate("text", x = 3, y = 4.5, label = "ANOVA, F(2)=49.16, p<0.001")+
scale_x_discrete(labels=c("I. setosa", "I. versicolor", "I. virginica"))
Comme pour les analyses univariées, cet atelier est seulement axé sur la création de graphiques. Pour une introduction aux analyses multivariées (ordinations, PERMANOVA, Procrustes), consulter l’atelier 6.
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-4
data(dune) # matrice de communautés
data(dune.env) # données environnementales
Ce jeu de données contient une matrice d’abondance (dune) d’espèces végétales avec leur classe de couverture pour 20 sites. Le dataframe (dune.env) contient 5 variables denvironnementales.
Transformation Hellinger
dune.hel <- decostand(dune, method="hellinger")
Faire une PCA (avec la fonction prcomp cette fois-ci, car l’objet qu’on va créer peut être “lu” plus tard par la fonction ggord)
pca.dune <- prcomp(dune.hel)
summary(pca.dune)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 0.4122 0.3333 0.2374 0.2097 0.20242 0.17376 0.14506
## Proportion of Variance 0.3057 0.1998 0.1013 0.0791 0.07371 0.05431 0.03785
## Cumulative Proportion 0.3057 0.5054 0.6068 0.6859 0.75960 0.81391 0.85177
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 0.13392 0.11745 0.10441 0.09995 0.09550 0.07985 0.07188
## Proportion of Variance 0.03226 0.02481 0.01961 0.01797 0.01641 0.01147 0.00929
## Cumulative Proportion 0.88403 0.90885 0.92846 0.94643 0.96283 0.97430 0.98360
## PC15 PC16 PC17 PC18 PC19 PC20
## Standard deviation 0.05424 0.05240 0.04063 0.03373 0.02536 1.816e-17
## Proportion of Variance 0.00529 0.00494 0.00297 0.00205 0.00116 0.000e+00
## Cumulative Proportion 0.98889 0.99383 0.99680 0.99884 1.00000 1.000e+00
Dans l’atelier 6, nous avons fait nos PCA avec la fonction rda. Lorsqu’on compare le sommaire des deux analyses, on voit que les deux fonctions donnent les mêmes résultats.
pca2.dune <- rda(dune.hel)
summary(pca2.dune)
##
## Call:
## rda(X = dune.hel)
##
## Partitioning of variance:
## Inertia Proportion
## Total 0.5559 1
## Unconstrained 0.5559 1
##
## Eigenvalues, and their contribution to the variance
##
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Eigenvalue 0.1699 0.1111 0.05634 0.04397 0.04097 0.03019 0.02104
## Proportion Explained 0.3057 0.1998 0.10135 0.07910 0.07371 0.05431 0.03785
## Cumulative Proportion 0.3057 0.5054 0.60679 0.68590 0.75960 0.81391 0.85177
## PC8 PC9 PC10 PC11 PC12 PC13
## Eigenvalue 0.01794 0.01379 0.01090 0.009991 0.00912 0.006376
## Proportion Explained 0.03226 0.02481 0.01961 0.017972 0.01641 0.011469
## Cumulative Proportion 0.88403 0.90885 0.92846 0.946428 0.96283 0.974303
## PC14 PC15 PC16 PC17 PC18 PC19
## Eigenvalue 0.005166 0.002942 0.002745 0.001651 0.001138 0.0006432
## Proportion Explained 0.009293 0.005292 0.004939 0.002969 0.002047 0.0011570
## Cumulative Proportion 0.983597 0.988888 0.993827 0.996796 0.998843 1.0000000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 1.80276
##
##
## Species scores
##
## PC1 PC2 PC3 PC4 PC5 PC6
## Achimill -0.2241355 0.057048 -0.068674 -0.173281 -0.0498964 -0.0006022
## Agrostol 0.4142446 -0.180068 0.020078 -0.020756 0.0029054 0.0193453
## Airaprae -0.0350694 0.154644 0.116317 -0.039999 -0.1384374 -0.0376930
## Alopgeni 0.1386550 -0.318444 0.186419 0.001329 0.0030171 -0.0215368
## Anthodor -0.1962910 0.238579 0.095768 -0.179434 -0.0380470 -0.0471885
## Bellpere -0.1148332 -0.054380 -0.057640 0.064957 0.0010453 0.0896219
## Bromhord -0.1410948 -0.053134 -0.049586 -0.057365 0.0210053 0.0938974
## Chenalbu 0.0121539 -0.025919 0.036102 -0.015353 -0.0038179 0.0173353
## Cirsarve -0.0026599 -0.031975 0.004872 0.021813 -0.0182940 0.0064877
## Comapalu 0.1100812 0.048635 -0.061981 -0.028257 0.0102828 0.1064664
## Eleopalu 0.3796078 0.069099 -0.185862 -0.064110 0.0144632 0.0035959
## Elymrepe -0.1456860 -0.248314 -0.094249 0.010195 -0.1482516 -0.0563573
## Empenigr 0.0004895 0.057213 0.063829 0.026314 -0.0393467 -0.0066958
## Hyporadi -0.0559119 0.197829 0.139068 0.036814 -0.1427093 -0.0326104
## Juncarti 0.2478729 -0.012497 -0.092252 0.005525 0.0265396 -0.1569375
## Juncbufo 0.0190348 -0.121457 0.179346 -0.051491 0.0931930 -0.0211675
## Lolipere -0.3366259 -0.173707 -0.206664 0.154192 0.0092155 -0.0325721
## Planlanc -0.2547353 0.197554 -0.034482 -0.039980 0.1469090 -0.0256698
## Poaprat -0.2808972 -0.142814 -0.082463 0.084918 -0.0502270 -0.0278725
## Poatriv -0.1317430 -0.363370 0.059442 -0.137436 0.0581366 -0.0136469
## Ranuflam 0.2783344 0.024399 -0.081013 -0.052662 0.0009058 0.0274683
## Rumeacet -0.1079283 -0.024094 0.054027 -0.122507 0.2108960 -0.0801798
## Sagiproc 0.0409392 -0.082141 0.246528 0.127747 -0.0148906 -0.0086761
## Salirepe 0.0651846 0.156446 0.007210 0.136628 -0.0166452 -0.0581183
## Scorautu -0.0562239 0.158686 0.094988 0.095244 0.0265794 0.1029567
## Trifprat -0.1001443 0.028383 -0.024352 -0.095019 0.1389927 -0.0417313
## Trifrepe -0.0471194 -0.006999 0.052812 0.045812 0.1395303 0.2537096
## Vicilath -0.0541878 0.052774 -0.023104 0.114929 0.0381791 0.0343033
## Bracruta 0.0852752 0.107955 -0.003783 0.181534 0.2198174 -0.1398609
## Callcusp 0.2045115 0.057120 -0.104072 -0.065216 -0.0146101 0.0579472
##
##
## Site scores (weighted sums of species scores)
##
## PC1 PC2 PC3 PC4 PC5 PC6
## 1 -0.423423 -0.42754 -6.533e-01 0.01583 -0.761705 -0.463469
## 2 -0.356333 -0.36816 -1.991e-01 -0.08370 -0.355682 0.630146
## 3 -0.080900 -0.55215 -5.753e-02 0.27275 -0.217514 0.016724
## 4 -0.041005 -0.49292 7.510e-02 0.33626 -0.282018 0.100013
## 5 -0.458003 0.05220 -1.432e-01 -0.54110 0.351891 -0.117334
## 6 -0.389342 0.20337 -1.022e-01 -0.30549 0.758789 -0.245696
## 7 -0.451813 0.06864 -6.835e-02 -0.41821 0.585536 -0.138732
## 8 0.295113 -0.29375 -6.362e-02 0.19828 0.010903 -0.185034
## 9 0.035973 -0.49264 2.714e-01 0.07945 0.086166 -0.419257
## 10 -0.453094 0.13639 -2.391e-01 -0.12201 0.144355 0.462946
## 11 -0.273214 0.29633 3.931e-05 0.87663 0.126762 0.097096
## 12 0.246622 -0.33173 9.205e-01 -0.03527 0.493526 -0.017653
## 13 0.226908 -0.48390 6.740e-01 -0.28663 -0.071279 0.323642
## 14 0.570818 0.25288 -3.165e-01 -0.32869 -0.085439 1.196310
## 15 0.654415 0.28845 -3.733e-01 0.01035 0.196968 0.002255
## 16 0.720030 -0.10689 -2.727e-01 -0.35992 0.007984 -0.504518
## 17 -0.317465 0.75272 3.395e-01 -0.64284 -0.803236 -0.262492
## 18 -0.201131 0.39818 -2.007e-01 0.89875 0.365682 0.086314
## 19 0.006263 0.73205 8.167e-01 0.33669 -0.503443 -0.085674
## 20 0.689579 0.36848 -4.077e-01 0.09888 -0.048247 -0.475587
PERMANOVA. On va tester si la composition végétale diffère entre différents type d’aménagement (facteur; variable quantitative) des sites.
dune.dist <- vegdist(decostand(dune, method='hellinger'), method='euclid')
disper.dune <- betadisper(dune.dist, dune.env$Management)
anova(disper.dune)
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.17812 0.059373 1.5729 0.2349
## Residuals 16 0.60397 0.037748
permanova.dune <- adonis2(dune.dist ~ Management, data = dune.env)
permanova.dune
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dune.dist ~ Management, data = dune.env)
## Df SumOfSqs R2 F Pr(>F)
## Management 3 3.1672 0.29986 2.2842 0.005 **
## Residual 16 7.3950 0.70014
## Total 19 10.5621 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Le test montre que le type d’aménagement a un effet statistiquement significatif sur la composition en espèces végétales des sites.
On va créer un graphique pour présenter la PCA avec les différents types d’aménagement.
Extraire les scores des sites et des espèces, à partir de l’objet créer par prcomp. Cela va nous permettre de mettre éventuellement seulement les espèces désirées sur le graphique.
site.scores<-pca.dune$x
site.scores.df <- data.frame(site.scores)
sp.scores<-pca.dune$rotation
sp.scores.df <- data.frame(sp.scores)
Ouverture des librairies nécessaires à la création du graphique. Si problème avec installation de la librairie ggord : https://rdrr.io/github/fawda123/ggord/
library(ggord)
library(ggplot2)
library(ggrepel) # pour éviter overlap du texte
On va d’abord commencer le graphique avec la fonctions ggord du package du même nom.
Le graphique créé par la fonction ggord est un objet ggplot. Cet objet (ici appelé pca.plot) peut être personnalisé à même la fonction ggord, et peut aussi servir de base sur laquelle ajouter des couches avec ggplot.
pca.plot <- ggord(pca.dune, xlims=c(-1,1), ylims=c(-1,1)) # xlims et ylims définissent les limites dans l'affichage des axe (les valeurs entre parenthèses correspondent à l'échelle numérique des axes). On ajuste au fur et à mesure, au besoin.
pca.plot # À la fin de chaque chunk, l'objet graphique est appelé pour qu'on visualise le rendu à chaque fois.
On remarque que la fonction
ggord inclut par défaut la proportion de la variance expliquée par les axes. C’est super!
Par défaut, la fonction ggord va mettre en graphique les deux premiers axes de l’ordination. Mais si on s’intéresse à différents axes, on peut spécifier lesquels mettre en graphique.
pca.plot <- ggord(pca.dune, xlims=c(-1,1), ylims=c(-1,1),
axes=c(1,3)) # Ici, on met en graphique les axes 1 et 3.
pca.plot
On revient aux deux premiers axes. Faisons quelques petits changements. Les changements sont marqués d’un “#” à côté de la ligne.
pca.plot <- ggord(pca.dune, xlims=c(-1.2,1.4), ylims=c(-0.75,0.9), # ajustement des limites
axes=c(1,2), # retour aux axes 1 et 2
size=3, # diminution de la taille des points (sites)
arrow=NULL,# enlever les flèches (qui viennent par défaut) pour les espèces
labcol="forestgreen") # changer la couleur du texte pour les espèces
pca.plot
Faisons d’autres changements!
pca.plot <- ggord(pca.dune, xlims=c(-1.2,1.4), ylims=c(-0.75,0.9),
axes=c(1,2),
size=3,
obslab=TRUE, # pour passer de points pour les sites (par défaut) à texte (nom des sites)
arrow=NULL,
labcol="forestgreen") # changer la couleur du texte pour les espèces
pca.plot
Il y a beaucoup d’espèces affichées au milieu et c’est illisible. En plus, les espèces au milieu ne contribuent pas beaucoup à la distinction compositionnelle des sites. On va donc garder seulement les quelques espèces dont les scores ont la valeur absolue la plus élevée, pour chacun des deux axes, qu’on met en graphique (ici PC1 et PC2).
D’abord, on va créer un nouveau dataframe contenant seulement ces espèces, en utilisant quelques fonctions de la librairie dplyr.
A <- top_n(sp.scores.df, 3, PC1) # sélectionne les 3 espèces dont le score est le plus élevé sur l'axe 1.
B <- top_n(sp.scores.df, -3, PC1) # sélectionne les 3 espèces dont le score est le plus bas sur l'axe 1.
C <- top_n(sp.scores.df, 3, PC2) # sélectionne les 3 espèces dont le score est le plus élevé sur l'axe 2.
D <- top_n(sp.scores.df, -3, PC2) # sélectionne les 3 espèces dont le score est le plus bas sur l'axe 2.
sp.scores.skim <- bind_rows(A,B,C,D) # merge les sub dataframe A, B, C et D.
sp.scores.skim <- distinct(sp.scores.skim) # on ne garde que les rangées (espèces) qui sont uniques (pour éviter qu'une même espèce s'y retrouve en copies)
Mise en graphique de ces espèces seulement
pca.plot <- ggord(pca.dune, xlims=c(-1.2,1.4), ylims=c(-0.75,0.9),
axes=c(1,2),
size=3,
obslab=TRUE,
arrow=NULL,
txt=NULL, # enlever le texte (par défaut) pour les espèces
addpts=sp.scores.skim, # Ajouter les espèces aux plus hauts et faibles scores
addsize=3, # taille des points des espèces
addcol="forestgreen") # couleur des espèces
# on enlève l'argument labcol
pca.plot
pca.plot <- ggord(pca.dune, xlims=c(-1.2,1.4), ylims=c(-0.75,0.9),
axes=c(1,2),
size=3,
obslab=TRUE,
arrow=NULL,
txt=NULL,
addpts=sp.scores.skim,
addsize=4, # augmenter la taille du texte des espèces
ptslab=TRUE, # changer de points à texte pour les espèces
addcol="forestgreen")
pca.plot
On va ajouter les résultats d’un envfit. D’abord, on fait le test.
(pca.envfit <- envfit(pca.dune, dune.env))
##
## ***VECTORS
##
## PC1 PC2 r2 Pr(>r)
## A1 0.97222 0.23409 0.3622 0.008 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## PC1 PC2
## Moisture1 -0.3635 0.0257
## Moisture2 -0.2224 -0.0314
## Moisture4 0.1408 -0.3321
## Moisture5 0.4504 0.0872
## ManagementBF -0.3597 0.0173
## ManagementHF -0.1930 -0.0745
## ManagementNM 0.2330 0.3751
## ManagementSF 0.1077 -0.3217
## UseHayfield -0.0994 0.2242
## UseHaypastu -0.0203 -0.2180
## UsePasture 0.1716 0.0350
## Manure0 0.2330 0.3751
## Manure1 -0.2294 -0.0161
## Manure2 -0.2385 -0.0895
## Manure3 0.1969 -0.1644
## Manure4 -0.1812 -0.3955
##
## Goodness of fit:
## r2 Pr(>r)
## Moisture 0.5366 0.001 ***
## Management 0.4614 0.001 ***
## Use 0.1794 0.131
## Manure 0.4531 0.010 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
Pour incorporer les vecteurs résultant de ce test, il faut d’abord les extraire et les mettre dans un dataframe.
vect <- data.frame(pca.envfit[["vectors"]][["arrows"]])
vect$variable.env<-rownames(vect)
pca.plot <- ggord(pca.dune, xlims=c(-1.2,1.4), ylims=c(-0.75,0.9),
axes=c(1,2),
size=3,
obslab=TRUE,
arrow=NULL,
txt=NULL,
addpts=sp.scores.skim,
addsize=4,
ptslab=TRUE,
addcol="forestgreen")
pca.plot +
coord_fixed() + # ajout des vecteurs de envfit
geom_segment(data = vect,
aes(x = 0, xend = PC1, y = 0, yend = PC2),
arrow = arrow(length = unit(0.5, "cm")), colour = "darkgrey") +
geom_text_repel(data = vect, aes(x = PC1, y = PC2, label = variable.env, size = 3), show.legend = FALSE)
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
Maintenant, on va ajouter des ellipses autour des sites appartenant à des groupes similaires. On va aussi revenir à des symboles au lieu du nom des sites et changer la couleur du texte des espèces.
pca.plot <- ggord(pca.dune, xlims=c(-1.2,1.4), ylims=c(-0.75,0.9),
axes=c(1,2),
grp_in=dune.env$Management, # ajout des ellipses en fonction du type de Management
grp_title = "Management",
ellipse_pro=0.95, # valeur de confiance pour les ellipses
alpha_el=0.3, # transparence des ellipses
size=3,
obslab=FALSE, # changement de texte à points pour sites
arrow=NULL,
txt=NULL,
addpts=sp.scores.skim,
addsize=4,
ptslab=TRUE,
addcol="black", # couleur des espèces
alpha = 1)
pca.plot +
scale_shape_manual('Groups', values = c(15,16,17,18)) + # Symboles pour les sites en fonction de leur groupe. Assurer vous d'avoir un nombre de symboles correspondant au nombre de groupes.
coord_fixed() +
geom_segment(data = vect,
aes(x = 0, xend = PC1, y = 0, yend = PC2),
arrow = arrow(length = unit(0.5, "cm")), colour = "darkgrey") +
geom_text_repel(data = vect, aes(x = PC1, y = PC2, label = variable.env, size = 3), show.legend = FALSE)
## Scale for shape is already present.
## Adding another scale for shape, which will replace the existing scale.
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
Voici une légende indiquant quels symboles sont associées aux différentes valeurs numériques de la couche scale_shape_manual.
Changer les couleurs des groupes On va extraire les codes des couleurs d’une palette de RColorBrewer avec la fonction brewer.pal
library(RColorBrewer)
pca.colors <- brewer.pal(n = 4, name = "Dark2") # Hexadecimal color specification
On va incorporer ces codes de couleurs dans un argument. On déplace la légende vers le haut.
pca.plot <- ggord(pca.dune, xlims=c(-1.2,1.4), ylims=c(-0.75,0.9),
axes=c(1,2),
grp_in=dune.env$Management,
grp_title = "Management",
cols=pca.colors, # ajout des couleurs choisies
ellipse_pro=0.95,
alpha_el=0.3,
size=3,
obslab=FALSE,
arrow=NULL,
txt=NULL,
addpts=sp.scores.skim,
addsize=4,
ptslab=TRUE,
addcol="black",
alpha = 1)
pca.plot <- pca.plot + # ici j'enregistre l'ensemble des arguments précédents et suivants dans mon objet pca.plot.
scale_shape_manual('Groups', values = c(15,16,17,18)) +
coord_fixed() +
geom_segment(data = vect,
aes(x = 0, xend = PC1, y = 0, yend = PC2),
arrow = arrow(length = unit(0.5, "cm")), colour = "darkgrey") +
geom_text_repel(data = vect, aes(x = PC1, y = PC2, label = variable.env, size = 3), show.legend = FALSE)+
theme(legend.position = 'top') # mettre la légende en haut
## Scale for shape is already present.
## Adding another scale for shape, which will replace the existing scale.
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
pca.plot
Créer un thème custom pour ggplot.
super_theme <- theme(
panel.background = element_blank(),
#panel.border = element_blank(),
panel.grid = element_blank(),
axis.line = element_line("gray25"),
text = element_text(size = 12),
axis.text = element_text(size = 10, colour = "gray25"),
axis.title = element_text(size = 14, colour = "gray25"),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14),
legend.key = element_blank())
Ajouter ce thème au graphique.
pca.plot + super_theme
On va rajouter du texte!
pca.plot + annotate("text", x = 1.1, y = 0.7, label = "PERMANOVA\nR2=30%\n p=0.007") +
super_theme
Satisfait? On sauvegarde l’image en png. Par défaut, ggplot enregistre le dernier graphique créé. Ici on copie-colle le code du graphique ci haut, puis sauvegarde ce graphique avec la fonction ggsave dans un dossier qu’on a créé (figures/).
pca.plot + annotate("text", x = 1.1, y = 0.7, label = "PERMANOVA\nR2=30%\n p=0.007") +
super_theme
ggsave("pca.png", path="./figures/")
## Saving 7 x 5 in image
Cet atelier n’est pas fait pour vous montrer à utiliser une esthétique particulière, mais bien pour vous familiariser avec les options (infinies…) permettant de personnaliser un graphique d’ordination. Vous allez voir des figures d’ordination dans la littérature. Essayez de reproduire l’esthétique qui vous semble la plus appropriée… ou, tout simplement, amusez-vous les artistes!!
Faire une CA avec la fonction cca.
ca.dune <- cca(dune)
ca.summary <- summary(ca.dune)
summary(ca.dune)
##
## Call:
## cca(X = dune)
##
## Partitioning of scaled Chi-square:
## Inertia Proportion
## Total 2.115 1
## Unconstrained 2.115 1
##
## Eigenvalues, and their contribution to the scaled Chi-square
##
## Importance of components:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7
## Eigenvalue 0.5360 0.4001 0.2598 0.17598 0.14476 0.10791 0.09247
## Proportion Explained 0.2534 0.1892 0.1228 0.08319 0.06844 0.05102 0.04372
## Cumulative Proportion 0.2534 0.4426 0.5654 0.64858 0.71702 0.76804 0.81175
## CA8 CA9 CA10 CA11 CA12 CA13 CA14
## Eigenvalue 0.08091 0.07332 0.05630 0.04826 0.04125 0.03523 0.020529
## Proportion Explained 0.03825 0.03466 0.02661 0.02282 0.01950 0.01665 0.009705
## Cumulative Proportion 0.85000 0.88467 0.91128 0.93410 0.95360 0.97025 0.979955
## CA15 CA16 CA17 CA18 CA19
## Eigenvalue 0.014911 0.009074 0.007938 0.007002 0.003477
## Proportion Explained 0.007049 0.004290 0.003753 0.003310 0.001644
## Cumulative Proportion 0.987004 0.991293 0.995046 0.998356 1.000000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
##
##
## Species scores
##
## CA1 CA2 CA3 CA4 CA5 CA6
## Achimill -0.90859 0.08461 -0.58636 -0.008919 -0.660183 -0.18877
## Agrostol 0.93378 -0.20651 0.28165 0.024293 -0.139326 -0.02256
## Airaprae -1.00434 3.06749 1.33773 0.194305 -1.081813 -0.53699
## Alopgeni 0.40088 -0.61839 0.85013 0.346740 0.016509 0.10169
## Anthodor -0.96676 1.08361 -0.17188 0.459788 -0.607533 -0.30425
## Bellpere -0.50018 -0.35503 -0.15239 -0.704153 -0.058546 0.07308
## Bromhord -0.65762 -0.40634 -0.30685 -0.496751 -0.561358 0.07004
## Chenalbu 0.42445 -0.84402 1.59029 1.248755 -0.207480 0.87566
## Cirsarve -0.05647 -0.76398 0.91793 -1.175919 -0.384024 -0.13985
## Comapalu 1.91690 0.52150 -1.18215 -0.021738 -1.359988 1.31207
## Eleopalu 1.76383 0.34562 -0.57336 -0.002976 -0.332396 -0.14688
## Elymrepe -0.37074 -0.74148 0.26238 -0.566308 -0.270122 -0.72624
## Empenigr -0.69027 3.26420 1.95716 -0.176936 -0.073518 -0.16083
## Hyporadi -0.85408 2.52821 1.13951 -0.175115 -0.311874 0.11177
## Juncarti 1.27580 0.09963 -0.09320 0.005536 0.289410 -0.78247
## Juncbufo 0.08157 -0.68074 1.00545 1.078390 0.268360 0.24168
## Lolipere -0.50272 -0.35955 -0.21821 -0.474727 0.101494 -0.01594
## Planlanc -0.84058 0.24886 -0.78066 0.371149 0.271377 0.11989
## Poaprat -0.38919 -0.32999 -0.02015 -0.358371 0.079296 -0.05165
## Poatriv -0.18185 -0.53997 0.23388 0.178834 -0.155342 -0.07584
## Ranuflam 1.55886 0.30700 -0.29765 0.046974 -0.008747 -0.14744
## Rumeacet -0.65289 -0.25525 -0.59728 1.160164 0.255849 -0.32730
## Sagiproc 0.00364 0.01719 1.11570 0.066981 0.186654 0.32463
## Salirepe 0.61035 1.54868 0.04970 -0.607136 1.429729 -0.55183
## Scorautu -0.19566 0.38884 0.03975 -0.130392 0.141232 0.23717
## Trifprat -0.88116 -0.09792 -1.18172 1.282429 0.325706 -0.33388
## Trifrepe -0.07666 -0.02032 -0.20594 0.026462 -0.186748 0.53957
## Vicilath -0.61893 0.37140 -0.46057 -1.000375 1.162652 1.44971
## Bracruta 0.18222 0.26477 -0.16606 0.064009 0.576334 0.07741
## Callcusp 1.95199 0.56743 -0.85948 -0.098969 -0.556737 0.23282
##
##
## Site scores (weighted averages of species scores)
##
## CA1 CA2 CA3 CA4 CA5 CA6
## 1 -0.81167 -1.0826714 -0.14479 -2.10665 -0.39287 -1.83462
## 2 -0.63268 -0.6958357 -0.09708 -1.18695 -0.97686 0.06575
## 3 -0.10148 -0.9128732 0.68815 -0.68137 -0.08709 -0.28678
## 4 -0.05647 -0.7639784 0.91793 -1.17592 -0.38402 -0.13985
## 5 -0.95293 -0.1846015 -0.95609 0.86853 -0.34552 -0.98333
## 6 -0.85633 -0.0005408 -1.39735 1.59909 0.65494 -0.19386
## 7 -0.87149 -0.2547040 -0.86830 0.90468 0.17385 -0.03446
## 8 0.76268 -0.2968459 0.35648 -0.10772 0.17507 -0.36444
## 9 0.09693 -0.7864314 0.86492 0.40090 0.28704 -1.02783
## 10 -0.87885 -0.0353136 -0.82987 -0.68053 -0.75438 0.81070
## 11 -0.64223 0.4440332 -0.17371 -1.09684 1.37462 2.00626
## 12 0.28557 -0.6656161 1.64423 1.71496 0.65381 1.17376
## 13 0.42445 -0.8440195 1.59029 1.24876 -0.20748 0.87566
## 14 1.91996 0.5351062 -1.39863 -0.08575 -2.21317 2.43044
## 15 1.91384 0.5079036 -0.96567 0.04227 -0.50681 0.19370
## 16 2.00229 0.1090627 -0.33414 0.33760 -0.50097 -0.76159
## 17 -1.47545 2.7724102 0.40859 0.75117 -2.59425 -1.10122
## 18 -0.31241 0.6328355 -0.66501 -1.12728 2.65575 0.97565
## 19 -0.69027 3.2642026 1.95716 -0.17694 -0.07352 -0.16083
## 20 1.94438 1.0688809 -0.66595 -0.55317 1.59606 -1.70292
Extraction de la proportion de la variance expliquée par les deux premiers axes (pour s’y référer plus rapidement).
prop.expl.ca1 <- ca.summary$cont$importance[2, "CA1"]
prop.expl.ca2 <- ca.summary$cont$importance[2, "CA2"]
envfit
(ca.envfit <- envfit(ca.dune, dune.env))
##
## ***VECTORS
##
## CA1 CA2 r2 Pr(>r)
## A1 0.998160 0.060614 0.3104 0.044 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## CA1 CA2
## Moisture1 -0.7484 -0.1423
## Moisture2 -0.4652 -0.2156
## Moisture4 0.1827 -0.7315
## Moisture5 1.1143 0.5708
## ManagementBF -0.7258 -0.1413
## ManagementHF -0.3867 -0.2960
## ManagementNM 0.6517 1.4405
## ManagementSF 0.3376 -0.6761
## UseHayfield -0.2861 0.6488
## UseHaypastu -0.0735 -0.5602
## UsePasture 0.5163 0.0508
## Manure0 0.6517 1.4405
## Manure1 -0.4639 -0.1738
## Manure2 -0.5872 -0.3600
## Manure3 0.5187 -0.3172
## Manure4 -0.2059 -0.8775
##
## Goodness of fit:
## r2 Pr(>r)
## Moisture 0.4113 0.005 **
## Management 0.4441 0.001 ***
## Use 0.1845 0.076 .
## Manure 0.4552 0.006 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
Pour incorporer les vecteurs résultant de ce test, il faut d’abord les extraire et les mettre dans un dataframe.
vect.ca.env <- data.frame(ca.envfit[["vectors"]][["arrows"]])
vect.ca.env$variable.env<-rownames(vect.ca.env)
PERMANOVA. On va tester si la composition végétale diffère entre différents type d’aménagement (facteur; variable quantitative) des sites.
dune.dist.chi2 <- vegdist(dune, method='chisq')
disper.dune.chi2 <- betadisper(dune.dist.chi2, dune.env$Management)
anova(disper.dune.chi2)
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 3.0277 1.0092 5.9124 0.006493 **
## Residuals 16 2.7311 0.1707
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova.dune.ca <- adonis2(dune.dist.chi2 ~ Management, data = dune.env)
permanova.dune.ca
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dune.dist.chi2 ~ Management, data = dune.env)
## Df SumOfSqs R2 F Pr(>F)
## Management 3 12.399 0.25508 1.8263 0.006 **
## Residual 16 36.208 0.74492
## Total 19 48.606 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Ici, le résultat du test de dispersion indique que la différence dans la composition végétale entre les types d’aménagement qui est “perçue” par la PERMANOVA pourrait être du à la différence dans la variance interne en composition des groupes.
Extraire les scores de CA avec la fonction scores
site.scores.dune.ca <- scores(ca.dune, display="sites", choices=c(1,2)) # ici on garde que les scores des axes 1 et 2 mais on aurait pu extraire tous les scores aussi.
site.scores.dune.ca <- data.frame(site.scores.dune.ca)
Management <- dune.env$Management
site.scores.ca <- cbind(site.scores.dune.ca, Management)
sp.scores.dune.ca <- scores(ca.dune, display="species", choices=c(1,2)) # idem scores axes 1 et 2 seulement
sp.scores.dune.ca <- data.frame(sp.scores.dune.ca)
Comme on a fait acvec les scores des espèces de la PCA, on ne garde que les plus hautes valeurs absolues.
A.ca <- top_n(sp.scores.dune.ca, 3, CA1)
B.ca <- top_n(sp.scores.dune.ca, -3, CA1)
C.ca <- top_n(sp.scores.dune.ca, 3, CA2)
D.ca <- top_n(sp.scores.dune.ca, -3, CA2)
sp.scores.skim.ca <- bind_rows(A.ca,B.ca,C.ca,D.ca)
sp.scores.skim.ca <- distinct(sp.scores.skim.ca)
Graphique simple avec les fonctions plot et ordiplot.
plot(ca.dune)
ordiellipse(ca.dune, dune.env$Management, display = "sites", kind = "sd", label = T)
On commence le graphique ggplot. On peut faire un graphique qui ressemble à celui fait par plot en ajoutant “à la main” les sites et les espèces (en fonction de leurs scores extraits).
ca.plot <-
ggplot(data = site.scores.ca, aes(x = CA1, y = CA2)) +
expand_limits(x=c(-4,4), y=c(-1,3)) + # sert à modifier l'échelle numérique du graphique
geom_text(data = site.scores.dune.ca, aes(x = CA1, y = CA2, label= rownames(site.scores.dune.ca)), color = "black") + # ajout du texte des sites
geom_text(data = sp.scores.dune.ca, aes(x = CA1, y = CA2, label= rownames(sp.scores.dune.ca)), color = "red") + # ajout du texte des espèces
xlab("CA1 (25.34%)") + ylab("CA2 (18.92%)") # contrairement avec la fonction ggord, on écrit la variation expliqué par les axes soi même
ca.plot
On peut rajouter les vecteurs de note envfit.
ca.plot <-
ggplot(data = site.scores.ca, aes(x = CA1, y = CA2)) +
expand_limits(x=c(-2,2.5), y=c(-1,3)) + # on zoom in sur l'échelle
geom_text(data = site.scores.dune.ca, aes(x = CA1, y = CA2, label= rownames(site.scores.dune.ca)), color = "black") +
geom_text(data = sp.scores.dune.ca, aes(x = CA1, y = CA2, label= rownames(sp.scores.dune.ca)), color = "red") +
xlab("CA1 (25.34%)") + ylab("CA2 (18.92%)") +
geom_segment(data = vect.ca.env, # ajout des vecteurs envfit
aes(x = 0, xend = CA1, y = 0, yend = CA2),
arrow = arrow(length = unit(0.75, "cm")), color = "black") +
geom_text_repel(data = vect.ca.env, aes(x = CA1, y = CA2, label = variable.env, size = 3), color = "firebrick1", fontface = "bold", show.legend = FALSE)
ca.plot
On va mettre des points pour les sites et changer leurs couleurs et leurs symboles en fonction de leur type d’aménagement. On va aussi mettre en graphique uniquement les espèces aux plus grands scores.
ca.plot <-
ggplot(data = site.scores.ca, aes(x = CA1, y = CA2), color = Management) + # on ajoute ici arg "color"
expand_limits(x=c(-2,2.5), y=c(-1,3)) +
geom_point(aes(shape = Management, color = Management), size=3) + # on change de geom_text à geom_point pour les sites
xlab("CA1 (25.34%)") + ylab("CA2 (18.92%)") +
geom_text_repel(data = sp.scores.skim.ca, aes(x = CA1, y = CA2), label=rownames(sp.scores.skim.ca), color = "black") + # on garde geom_text pour les espèces mais on change de dataframe
geom_segment(data = vect.ca.env, # ajout des vecteurs envfit
aes(x = 0, xend = CA1, y = 0, yend = CA2),
arrow = arrow(length = unit(0.75, "cm")), color = "black") +
geom_text_repel(data = vect.ca.env, aes(x = CA1, y = CA2, label = variable.env, size = 3), color = "firebrick1", fontface = "bold", show.legend = FALSE) +
scale_shape_manual(values = c(15,16,17,18)) # on décide des symboles à utiliser (se référer à la charte plus haut)
ca.plot
On va ajouter des ellipses autour de nos sites en fonction de leur type d’aménagement. On va d’abord utiliser la fonction ggordiplot de la librairie ggordiplots pour extraire les données des ellipses. (on peut aussi extraire les données de hulls et spiders.) En même temps, ça nous donne un autre beau graphique, mais qui est plus complexe à modifier.
library(ggordiplots)
## Loading required package: glue
ca.ggordiplot <- gg_ordiplot(ca.dune, groups = dune.env$Management, kind = "sd", conf = 0.95, pt.size = 3) # on spécifie l'intervalle de confiance des ellipses
names(ca.ggordiplot)
## [1] "df_ord" "df_mean.ord" "df_ellipse" "df_hull" "df_spiders"
## [6] "plot"
ca.ellipses <- ca.ggordiplot$df_ellipse # extraction des valeurs des ellipses
On prend les valeurs des ellipses qu’on vient d’extraire pour mettre en graphique dans ggplot ces ellipses. Om change aussi les couleurs de nos groupes.
ca.plot <-
ggplot(data = site.scores.ca, aes(x = CA1, y = CA2), color = Management) +
expand_limits(x=c(-2,2.5), y=c(-1,3)) +
geom_point(aes(shape = Management, color = Management), size=3) +
xlab("CA1 (25.34%)") + ylab("CA2 (18.92%)") +
geom_text_repel(data = sp.scores.skim.ca, aes(x = CA1, y = CA2), label=rownames(sp.scores.skim.ca), color = "black") +
geom_segment(data = vect.ca.env, # ajout des vecteurs envfit
aes(x = 0, xend = CA1, y = 0, yend = CA2),
arrow = arrow(length = unit(0.75, "cm")), color = "black") +
geom_text_repel(data = vect.ca.env, aes(x = CA1, y = CA2, label = variable.env, size = 3), color = "firebrick1", fontface = "bold", show.legend = FALSE) +
scale_shape_manual(values = c(15,16,17,18)) +
geom_path(data = ca.ellipses, aes(x = x, y = y, color = Group), show.legend = FALSE) + # ajout des ellipses
scale_color_brewer(palette="Dark2") # changement des couleurs des ellipses et des sites
ca.plot
Enfin, on pourrait rajouter notre thème personnalisé, mais optons pour un existant comme “bw”.
ca.plot + theme_bw()
Comme mentionnée précédemment, vous pouvez modifier presque à l’infini votre graphique!
On sauvegarde notre figure de CA.
ca.plot + theme_bw()
ggsave("ca.png", path="./figures/")
## Saving 7 x 5 in image
Faire une PCoA avec la fonction prcomp.
dune.bray <- vegdist(dune, method='bray')
pcoa.dune <- cmdscale(dune.bray, k =(nrow(dune) - 1), eig = TRUE)
## Warning in cmdscale(dune.bray, k = (nrow(dune) - 1), eig = TRUE): only 14 of the
## first 19 eigenvalues are > 0
summary(pcoa.dune)
## Length Class Mode
## points 280 -none- numeric
## eig 20 -none- numeric
## x 0 -none- NULL
## ac 1 -none- numeric
## GOF 2 -none- numeric
Extraire les scores de la PCoA.
site.scores.pcoa <- scores(pcoa.dune) # site scores
sp.scores.pcoa <- wascores(pcoa.dune$points, dune) # species scores
site.scores.pcoa <- data.frame(site.scores.pcoa)
sp.scores.pcoa <- data.frame(sp.scores.pcoa)
Management <- dune.env$Management
site.scores.pcoa <- cbind(site.scores.pcoa, Management)
Comme on a fait acec les scores des espèces de la PCA et de la CA, on ne garde que les plus hautes valeurs absolues, mais cette fois ci, des axes 1 et 3 (pour le plaisir)!.
A.pcoa <- top_n(sp.scores.pcoa, 3, X1)
B.pcoa <- top_n(sp.scores.pcoa, -3, X1)
C.pcoa<- top_n(sp.scores.pcoa, 3, X3)
D.pcoa <- top_n(sp.scores.pcoa, -3, X3)
sp.scores.skim.pcoa <- bind_rows(A.pcoa,B.pcoa,C.pcoa,D.pcoa)
sp.scores.skim.pcoa <- distinct(sp.scores.skim.pcoa)
envfit SUR LES AXES 1 ET 3 (on n’oublie pas!)
(pcoa.envfit <- envfit(pcoa.dune, dune.env, choices=c(1,3)))
##
## ***VECTORS
##
## Dim1 Dim3 r2 Pr(>r)
## A1 0.9999700 0.0078085 0.379 0.021 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## Dim1 Dim3
## Moisture1 -0.2650 0.0607
## Moisture2 -0.1620 0.0000
## Moisture4 0.1065 -0.2369
## Moisture5 0.3271 0.0069
## ManagementBF -0.2694 0.0706
## ManagementHF -0.1350 -0.0394
## ManagementNM 0.1822 0.0370
## ManagementSF 0.0650 -0.0395
## UseHayfield -0.0608 -0.0240
## UseHaypastu -0.0227 -0.0210
## UsePasture 0.1213 0.0673
## Manure0 0.1822 0.0370
## Manure1 -0.1654 0.0175
## Manure2 -0.1648 -0.0944
## Manure3 0.1397 -0.0451
## Manure4 -0.1656 0.0944
##
## Goodness of fit:
## r2 Pr(>r)
## Moisture 0.6918 0.001 ***
## Management 0.2634 0.141
## Use 0.0614 0.640
## Manure 0.2892 0.203
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
Pour incorporer les vecteurs résultant de ce test, il faut d’abord les extraire et les mettre dans un dataframe.
vect.pcoa.env <- data.frame(pcoa.envfit[["vectors"]][["arrows"]])
vect.pcoa.env$variable.env <- rownames(vect.pcoa.env)
PERMANOVA. On va tester si la composition végétale diffère entre différents type d’aménagement (facteur; variable quantitative) des sites.
dune.dist.bray <- vegdist(dune, method='bray')
disper.dune.bray <- betadisper(dune.dist.bray, dune.env$Management)
anova(disper.dune.bray)
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.13831 0.046104 1.9506 0.1622
## Residuals 16 0.37816 0.023635
permanova.dune.pcoa <- adonis2(dune.dist.bray ~ Management, data = dune.env)
permanova.dune.pcoa
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dune.dist.bray ~ Management, data = dune.env)
## Df SumOfSqs R2 F Pr(>F)
## Management 3 1.4686 0.34161 2.7672 0.001 ***
## Residual 16 2.8304 0.65839
## Total 19 4.2990 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Extraction des ellipses de la PCoA (axes 1 et 3!)
library(ggordiplots)
pcoa.ggordiplot <- gg_ordiplot(pcoa.dune, choices = c(1,3), groups = dune.env$Management, kind = "sd", conf = 0.95, pt.size = 3)
names(pcoa.ggordiplot)
## [1] "df_ord" "df_mean.ord" "df_ellipse" "df_hull" "df_spiders"
## [6] "plot"
pcoa.ellipses <- pcoa.ggordiplot$df_ellipse
pcoa.plot <-
ggplot(data = site.scores.pcoa, aes(x = Dim1, y = Dim3), color = Management) +
expand_limits(x=c(-0.35,0.6), y=c(-0.35,0.35)) +
geom_point(aes(shape = Management, color = Management), size=3) +
xlab("PCoA1 (X%)") + ylab("PCoA3 (X%)") +
geom_text_repel(data = sp.scores.skim.pcoa, aes(x = X1, y = X3), label=rownames(sp.scores.skim.pcoa), color = "black") +
geom_segment(data = vect.pcoa.env,
aes(x = 0, xend = Dim1, y = 0, yend = Dim3),
arrow = arrow(length = unit(0.75, "cm")), color = "black") +
geom_text_repel(data = vect.pcoa.env, aes(x = Dim1, y = Dim3, label = variable.env, size = 3), color = "firebrick1", fontface = "bold", show.legend = FALSE) +
scale_shape_manual(values = c(15,16,17,18)) +
geom_path(data = pcoa.ellipses, aes(x = x, y = y, color = Group), show.legend = FALSE) +
scale_color_brewer(palette="Dark2")
pcoa.plot
ggsave("pcoa.png", path="./figures/")
## Saving 7 x 5 in image